f = @(u,t) -5*u;
u_true = @(t) exp(-5*t);

k = 0.5;
t0 = 0;
U0 = 1;

% Step 1: Compute U*
% Implicit midpoint (solve for U*)
U_star_fun = @(U_star) U_star - (U0 + k/4 * (f(U0,t0) + f(U_star,t0 + k/2)));
U_star = fsolve(U_star_fun, U0);  % initial guess = U0

% Step 2: Compute U^{n+1}
% Implicit BDF2 (solve for U_{n+1})
U1_fun = @(U1) U1 - (1/3)*(4*U_star - U0 + k*f(U1, t0 + k));
U1 = fsolve(U1_fun, U0);

% Exact solution for plotting
tt = linspace(t0, t0 + k, 100);
uu = u_true(tt);

% Plot


figure; hold on;
plot(tt, uu, 'k', 'LineWidth', 1.5);  % true solution
plot(t0, U0, 'ko', 'MarkerFaceColor', 'k'); % start
plot(t0 + k/2, U_star, 'bo', 'MarkerFaceColor', 'b'); % midpoint
plot(t0 + k, U1, 'ro', 'MarkerFaceColor', 'r'); % end
plot(t0+2/3*k, 4/3*U_star - 1/3*U0, 'ko', 'MarkerFaceColor', 'g'); % end

% Arrows
quiver(t0, U0, k, k*f(U0,t0), 0, 'b:', 'LineWidth', 1, 'MaxHeadSize', 0)
%text(t0+k/5 , U0+ k/3*f(U0,t0), 'l_1 : f(U^n, t_n)', 'FontSize', 10, 'Color', 'k');
quiver(t0+k/2, U_star, k/2, k/2*f(U_star, t0+k/2), 0, 'b-.', 'LineWidth', 1, 'MaxHeadSize', 0)
%text(t0+k*0.5 ,U_star+0.02 , 'l_2 : f(U^*, t_n+k/2)', 'FontSize', 10, 'Color', 'b');
% quiver(t0+k/2, U1-k/2*f(U1, t0+k), k/2, k/2*f(U1, t0+k), 0, 'b--', 'LineWidth', 1, 'MaxHeadSize', 0)
quiver(t0+2/3*k, 4/3*U_star - 1/3*U0, k/3, k*f(U1, t0 + k)/3, 0, 'b-', 'LineWidth', 1, 'MaxHeadSize', 0)
%text(t0+k*0.8 , U0+ k*0.8*f(U0,t0), 'l_3 : f(U^{n+1}, t_{n+1})', 'FontSize', 10, 'Color', 'r');
quiver(t0+k, u_true(t0+k), 0, U1 - u_true(t0+k), 0, 'r', 'LineWidth', 2, 'MaxHeadSize', 0)

legend('True solution','U^n', 'U^*', 'U^{n+1}','4/3*U^* - 1/3*U^n','l_1','l_2','l_3','Location', 'best')
title('Geometric interpretation of TR-BDF2')
xlabel('t'), ylabel('u')
grid on; axis tight;

